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We present the Bethe ansatz solution for the discrete time zero range and asymmetric exclusion 
processes with fully parallel dynamics. The model depends on two parameters: p, the probability 
of single particle hopping, and q, the deformation parameter, which in the general case, |g| < 1, is 
responsible for long range interaction between particles. The particular case q = corresponds to 
the Nagel-Schreckenberg traffic model with «max = 1- As a result, we obtain the largest eigenvalue 
of the equation for the generating function of the distance travelled by particles. For the case 
q = the result is obtained for arbitrary size of the lattice and number of particles. In the general 
case we study the model in the scaling limit and obtain the universal form specific for the Kardar- 
Parisi-Zhang universality class. We describe the phase transition occurring in the limit p — > 1 when 
<? < 0. 
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I. INTRODUCTION. 



One-dimensional models of stochastic processes attracted much attention 
last decade. Being related to several natural phenomena like the interface 
growth 0, the traffic flowj^l, and the self-organized criticality 3] , they ad- 
mit an exact calculation of many physical quantities, which can not be ob- 
tained with mean-field approach. Such models served as a testing ground 
for the description of many interesting effects specific for nonequilibrium 
systems like, boundary induced phase transitions shock waves 6J 

and condensation transition . 

The most prominent in the physical community one-dimensional 
stochastic model is the Asymmetric Simple Exclusion Process (ASEP). 
The interest in this model was inspired by its Bethe ansatz solution 0,13, 
which became the first direct exact calculation of the dynamical exponent 
of Kardar-Parisi-Zhang equation T|. Since then plenty of exact results 
were obtained in this direction using both the Bethe ansatz and the 
matrix product method [T]| . The advantage of the former is the possibility 
of consideration of time dependent quantities rather than only the station- 
ary ones. Among the results obtained for the ASEP with the help of the 
Bethe Ansatz there are the crossover scaling functions for the Kardar- 
Parisi-Zhang universality class the large deviation function of the 
Kardar-Parisi-Zhang-type interface ^^e time dependent correla- 
tion fimctions on the infinite and the periodic lattices [l^. 



2 



Though the Bethe ansatz solvabihty opens rich opportunities for ob- 
taining exact results, it imphes restrictive hmitations for the dynamical 
rules, such that the range of the Bethe ansatz solvable models is not too 
wide. Besides the ASEP, the examples of the particle hopping models, 
which admit the Bethe ansatz solution, are the asymmetric diffusion mod- 
els ^3 I generalized drop-push models 20^ 21] , and the asymmetric 
avalanche process (ASAP) |23|. Most of the models of stochastic pro- 
cesses solved by the Bethe ansatz are formulated in terms of continuous 
time dynamics or random sequential update, which allows one to use the 
analogy with the integrable quantum chains. The results for the processes 
with discrete time parallel update are rare [l^. On the other hand in 
conventional theory of quantum integrable systems the fundamental role is 
played by the two dimensional vertex models, all the quantum spin chains 
and continuous quantum models being their particular limiting cases p^ . 
In the theory of one-dimensional stochastic processes the same role could 
be played by the models with discrete time parallel update, i.e. so called 
stochastic cellular automata [2^. Having plenty of real applications, the 
cellular automata give splendid opportunity for doing large scale numeri- 
cal simulations. Thus, finding the Bethe ansatz solutions for discrete time 
models with fully parallel update would be of interest. The lack of such 
solutions owes probably to more complicated structure of the stationary 
state of models with the fully parallel update, which makes the application 
of the Bethe ansatz more subtle. 

Recently the Bethe ansatz was applied to solve the continuous time 
zero-range process (ZRP) j2Qj with the nonuniform stationary state p^ . 
The solution was based on the Bethe ansatz weighted with the stationary 
weights of corresponding configurations. The simple structure of the sta- 
tionary state of the ZRP allowed finding the one-parametric family of the 
hopping rates, which provides the Bethe ansatz integrability of the model. 
In the present article we use the same trick to solve much more general 
model of the ZRP with fully parallel update, which is defined by two- 
parametric family of hopping probabilities. As particular limiting cases of 
the parameters the model includes q-boson asymmetric diffusion mo del [I^. 
which in turn includes the drop-push |19| and phase model [2^ . and the 
ASAP (23. Simple mapping allows one to consider also the ASEP-like 
model that obeys the exclusion rule. In general it includes the long range 
interaction between particles, which makes the hopping probabilities de- 
pendent on the length of queue of particles next to the hopping one. In 
the simplest case, when the latter interaction is switched off, the model 
reduces to the Nagel-Schreckenberg trafhc model 2] with the maximal ve- 
locity Wmax = 1. In the article we find the eigenfunctions and eigenvalues 
of the equation for the generating function of the distance travelled by par- 
ticles and obtain the Bethe equations for both ASEP and ZRP cases. We 
analyze the solution of the Bethe equations corresponding to the largest 
eigenvalue of the equation for the generating function of the distance trav- 
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elled by particles. The analysis shows that a particular limit of parameters 

of the model reveals the second order phase transition. When the den- 
sity of particles reaches the critical value, the new phase emerges, which 
changes the dependence of the average flow of particles on the density of 
particles. The phase transition turns out to be intimately related to the 
intermittent-to-continuous flow transition in the ASAP, and the jamming 
transition in the traffic models. We analyze this transition in detail. 

The article is organized as follows. In the section II we formulate the 
model and discuss the basic results of the article. In the section III we 
give the Bcthe ansatz solution of the equation for the generating function 
of the moments of the distance travelled by particles. In the section IV 
we study the long time behavior of this equation and obtain its largest 
eigenvalue. The special form of the Bethe equations in the particular case, 
which corresponds to the Nagel-Schreckenberg model with Umax = 1, allows 
one to proceed with the calculations for arbitrary size of the lattice and 
number of particles. In the other cases we obtain the results in the scaling 
limit. In the section V we discuss the particular limit of the model, which 
exhibits the phase transition. In the section VI we discuss the calculation of 
the stationary correlation functions with the partition function formalism. 
A short summary is given in the section VII. 

II. MODEL DEFINITION AND MAIN RESULTS. 

A. Zero range process. 

The system under consideration can be most naturally formulated in 
terms of the discrete time ZRP. Let us consider the one-dimensional pe- 
riodic lattice consisting of N sites with M particles on it. Every site can 
hold any number of particles. A particle from an occupied site jumps to 
the next site forward with probability p{n), which depends only on the 
occupation number n of a site of departure. The system evolves step by 
step in the discrete time t, all sites being updated simultaneously at every 
step. Thus, as a result of the update the conflguration C, deflned by the 
set of occupation numbers 

C = {m, . . . ,njv} , (1) 

changes as follows 

{m, . . . , njv} ^ {ni - k2 + ki, . . . ,nN - ki + fcjv} • (2) 

Here the variable ki , taking values 1 or 0, denotes the number of particles 
arriving at the site i. According to these dynamical rules, all ki-s asso- 
ciated to the same time step are independent random variables with the 
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distribution, which depends on the occupation number of the site i — 1 

P{h = l)=p ,P{h = 0) = l-p . (3) 

The probability Pt (C) for the system to be in a configuration C at time 
step t obeys the Markov equation 

P*+i(C) = 5]T(C,C")P*(C'), (4) 

{C} 

where T (C, C") is the probabihty of the transition from C" to C. This 
equation is known to have a unique stationary solution which belongs 
to the class of the so called product measures, i.e. the probability of a 
configuration is given by the product of one-site factors 

1 ^ 

where the one-site factor / (n) is defined as follows 
/(O) = 
fin) = 

and Z {N, M) is the normalization constant. The properties of this sta- 
tionary measure have been intensively investigated particularly due to the 
possibility of condensation transition in such systems ^ . The question we 
are going to address below is: which dynamics leading to such a stationary 
measure admits the exact solution? 



B. Associated exclusion process. 

A simple mapping allows one to convert ZRP into ASEP. Given a ZRP 
configuration, we represent every site occupied with n particles as a se- 
quence of n sites occupied each with one particle plus one empty site ahead 
(see Fig. P). Thus, we add M extra sites to the lattice, so that its size 
becomes L = M + N. Obviously the exclusion constraint is imposed now. 
The evolution started with such a configuration is completely defined by 
the above ZRP dynamics. Specifically, any particle, which belongs to a 
cluster of n particles and has empty site ahead (i.e. the first particle of 
the cluster), jumps forward with the probability p(n) or stays with the 
probability {l — p{n)). Thus, we obtain the generalized ASEP with the 
asymmetric long range interaction, which we refer to as the ASEP associ- 
ated to given ZRP. Note that the above ZRP-ASEP relation applied on the 
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FIG. 1: ZRP-ASEP mapping. Each site occupied with n particles is mapped to 
the sequence of n sites occupied each by one particle with one empty site ahead. 

periodic lattice establishes correspondence between the sequences of con- 
figurations rather than between the configurations themselves. This fact 
can be illustrated with a simple example. Consider the following sequence 
of steps of the ZRP evolution. At the first step a particle jumps from an 
occupied site i to the next site i + 1, at the second from i + 1 to i + 2 and 
so on until a particle comes to the site i from i — 1 restoring the initial 
configuration. Reconstruction of the corresponding sequence of the ASEP 
configurations implies that when a particle joins a cluster of particles from 
the behind, the front particle of the same cluster takes the next step, which 
leads to the shift of the cluster one step backward. As a result we finally 
come to the configuration translated as a whole one step backward with 
respect to the initial one. This difference does not allow the ZRP and the 
ASEP to be considered as one model. There are, however, some quantities, 
which are not sensitive to translations^, and therefore arc identical for two 
models. For example the stationary measure of the associated ASEP is 
given by the expression (O , though / (n) corresponds now to a cluster of 
n particles with one empty site ahead. Apparently, the probability of any 
configuration does not change under the translation of the configuration as 
a whole. Below we consider the generating function of the total distance 
travelled by particles, Yj, in the infinite time limit, Hmt^oo In (e''^*) /t 
,which is also the same for the ZRP and the ASEP associated to it. The 
only thing we should keep in mind in such cases is that in ZRP the numbers 
and M denote the lattice size and the number of particles respectively, 
while in the associated ASEP N is the number of holes, M is the number 
of particles, and the size of the lattice is L = iV + M. Of course, time- 
or spatially- dependent quantities in general require separate consideration 
for each model. In this article we first consider ZRP dynamics and then 



^ The exact meaning of the term "insensitive to translations" will be clarified in the 
section III. 



point out the difference of solution for the of associated ASEP. 

In the end of the discussion of the ZRP-ASEP correspondence we should 
mention that another version of the ASEP can be obtained from the above 
one by the particle-hole transformation. Then the hopping probabilities 
will depend on the length of headway in front of the particle like for ex- 
ample the hopping rates in the bus route model l^^. Though the latter 
formulation seems more naturally related to the real traffic, we will use 
the former one to keep the direct connection with the ZRP. Of course the 
results for both versions can be easily related to each other. 

C. The results. 

To formulate basic result of the article we introduce the following gener- 
ating function 

oo 

y=o 

where Pt (C, Y) is the joint probability for the system to be in configuration 
C at time t, the total distance travelled by particles Yt being equal to Y. By 
definition the generating function, Ft (C), coincides with the probability 
Pt (C) of the configuration C in the particular case 7 = 0. The function 
Ft{C) obeys the evolution equation similar to I^J 

Ft+i{C) - e^^(^^^')r (C, C) Ft (C) . (7) 

{C} 

The term e'''^('^''^ ) accounts the increase of the total distance Yt, N (C, C) 
being the number of particles, which make a step during the transition from 
C to C . Below we show that the eigenfunctions of this equation, satisfying 
the following eigenfunction problem 

A (7) F,{C) = (C, C) F^ (C) , (8) 

{C} 

have the form of the Bethe ansatz weighted with the weights of stationary 
configurations provided that the hopping probabilities has the following 
form 

p{n)=px[n]^, n = l,2,3,..., (9) 
where [n]^ are the so called q- numbers 
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Thus, as usual, the Bethe ansatz to be applicable, the mfinite in general set 
of toppling probabilities should be reduced to two parametric family, the 
parameters being p and q. Formally, there are no any further limitations 
on the parameters . However, once we want p (n) to be probabilities, they 
should be less than or equal to one and nonnegative. It is obviously correct 
for any n if 

0<P<l-q and l^l < 1. (11) 

If g > 1, and an integer n* > 1 exists, such that 1/ [n* + 1]^ < p < 1/ [n*] 
we can still formulate the model with a finite number of particles M < n . 
The only way to consider an arbitrary number of particles with g > 1 is 
to consider the limit p — > +0. Particularly, if we put p — St, where St 
is infinitesimally small, the model turns into the continuous time version 
of ZRP (or associated ASEP) known also as q-boson asymmetric diffusion 
model |T7|. In this model the Poisson rate u{n) of the hopping of a particle 
from a site is given by u{n) = [n]q, q taking on values in the range q G 
(—1, oo). Such continuous time limit seems to be the only possibility for q 
to exceed 1. This model was studied in |0|,[23|. The case {p > 1 — q, 
\q\ < 1) also implies p (n) > 1 for some finite n > n* and thus do not allow 
consideration of arbitrary M, which is of practical interest. Thus, below 
we concentrate on the domain (|11|) . 

The other particular limits of the model under consideration, which re- 
produce the models studied before are as follows: 

- The Nagel-Schreckenberg traffic model with Umax = 1 corresponds to 
the ASEP with q = 0, when the probabilities p{n) are independent of 
n, p (n) = p. This case will be considered separately because, due to 
the special factorized form of the Bethe equations, it can be treated for 
arbitrary finite TV and M. 

— The asymmetric avalanche process corresponds to the limit 

(l-p)-O. (12) 

Originally the ASAP has been formulated as follows. In a stable state, 
M particles are located on a ring of sites, each site being occupied at 
most by one particle. At any moment of continuous time each particle can 
jump one step forward with Poisson rate 1. If a site i contains more than 
one particle rii > 1, it becomes unstable, and must relax immediately by 
spilling forward either n particles with probability ^„ or n — 1 particles 
with probability 1 — /i„ . The relaxation stops when all sites become stable 
again with < 1 for any i. The period of subsequent relaxation events is 
called avalanche. The avalanche is implied to be infinitely fast with respect 
to the continuous time. The toppling probabilities, which ensure the exact 
solvability of the model can also be written in terms of g— numbers 



(13) 
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Below we show that the ASAP can be obtained as the continuous time 
hmit of the ZRP observed from the moving reference frame, one step of the 
discrete time t being associated with the infinitesimal time interval (1—p). 

Let us now turn to the physical results we can extract from the Bethe 
ansatz solution. To this end, we note that the large time behavior of the 
generating function of the distance travelled by particles (e"*"^') is defined 
by the largest eigenvalue Aq (7) of the equation (|SJ| 

(e^^')=5]F,(C)~(Ao(7))*, t^oo. 

{C} 

Respectively, the logarithm of the largest eigenvalue Aq (7) is the generating 
function of the cumulants of Yt in the infinite time limit 



hm ^^"^^ - 9"lnAo(7) 



t^oo t 97" 



7=0 



Here the angle brackets denote the cumulants rather than moments of the 
particle current distribution, which is emphasized with the subscript c. 
Below we give the results for Aq (7) given in terms of two parameters, one 
of which is q introduced above, and the other is A defined as follows 

A = , (14) 

which captures all the dependence on the above p. 

* For q = 0, which corresponds to the Nagel-Schreckenberg traffic model 
with unit maximal velocity, fmax = 1, we obtain 

, , , . .^B"fLn-2\ ^ f 1- Mn,l~nN ,\ 

lnAo(7) = -A}_^— ( J 2i^i ( o_„r ;-A),(15) 



n 

71 = 1 



\Nn -ij \ 2-nL 



' M ^ n \ Nn J \ I - nL 

where Aq (7) is defined parametrically, both Aq (7) and 7 being power series 

in the variable B, (^) — is a binomial coefficient, 2F1 ^ '^^^ 

is the Gauss hypergeometric function, and L = N + M . 

* For an arbitrary q we obtain the results in the scaling limit N — > 00, 
M 00, M/N = p = const, ^N^/'^ = ccmst 

In Ao(7) = N(j)^ + kiN-^l'^G{k2N^''^i), (17) 



where G{x) is a universal scaling function defined parametrically, both G{x) 
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and X being the functions of the same parameter >c 

oo g 

S — 1 

OO g 

- = -E^- (19) 

This function appeared before in the studies of the ASEP ASAP [s^ 
and polymers in random media [s^ beheved to belong to the KPZ univer- 
sality class. The parameters 0, ki, k2 are the model dependent constants, 
which depend on two parameters A and q and the density of particles p. 
An explicit form of these constants will be derived in the following sections. 

The results in the limit p ^ \,q < Q are of special interest. Below 
we show that in this limit the model exhibits a phase transition, which is 
closely related to the intermittent-continuous flow transition in the ASAP 
and the jamming transition in traffic models. When the density of particles 
is less than critical density 

Pc-l/(l-9), (20) 

almost all particles synchronously jump forward with rare exceptions, 
which happen with probability of order of (1 — p). Thus the average flow of 
particles is equal to the density of particles with small correction of order 
of 

0(p<Pc)-P (21) 
However, when the density approaches the average flow (j) stops growing. 

(p > Pc) - Pc (22) 

It turns out that a fraction of particles gets stuck at a small fraction of sites 
instead of being involved into the particle flow. As a result the density of 
the particles involved into the flow is kept equal to p^. This phenomenon 
can be treated as a phase separation. The new phase consisting of immobile 
particles emerges at the critical point. Any increase of the density above 
the critical point leads to the growth of this phase, while the density of 
mobile particles remains unchanged. In terms of the associated ASEP 
or the traffic models one can describe the situation as an emergence of a 
small number of traffic jams with average length diverging, when {\ — p) 
decreases. Such behaviour leads to singularities in the constants (/>, A:i,A:2i 
Eq. H17|l , which define the large deviations of the particle current from the 
average. Below we analyze them in detail. 

We also study the stationary measure of the model with the canonical 
partition function formalism. The calculation of the stationary correlation 
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functions are discussed. We observe the relations between the Bethe ansatz 
solution in the thermodynamic limit and the thermodynamic properties 
of the stationary state. We make a conjecture that the relations found 
between the parameters characterizing the large deviations of the particle 
current and the parameters of the stationary state hold for the general ZRP 
belonging to the KPZ universality class. Finally we analyze the change of 
the behaviour of the occupation number distribution at the critical point in 
the limit p ^ 1. In addition the evaluation of Z (N, M) reveals interesting 
relation of the model with the theory of g-series , being based on the 
use of the most famous relation of this theory, q-binomial theorem. 

III. BETHE ANSATZ SOLUTION. 

Let us consider the equation for the generating function Ft (C) for the 
ZRP, 

N 

Ft+i{ni,...,nN) = ]~[ (e> {uj - h + kj+i))^'*^ 

{k,}^=l 

X + (23) 

xFt (ni - fci + fc2, . . . ,njv - fcjv + fci) , 

where we imply periodic boundary conditions + 1 = 1 . The summation 
is taken over all possible values of ki, . . . ,kN, which are the numbers of 
particles arriving at sites 1, . . . , N respectively. They take on the values 
of either or 1 if 7^ and of only otherwise. 

Before using the Bethe ansatz for the solution of this equation we should 
make the following remark. Most of the models studied by the coordi- 
nate Bethe ansatz have a common property. That is, a system evolves to 
the stationary state, where all the particle configurations occur with the 
same probability. This property can be easily understood from the struc- 
ture of the Bethe eigenfunction. Indeed, the stationary state is given by 
the groundstate of the evolution operator, which is the eigenfunction with 
zero eigenvalue and momentum. Such Bethe function does not depend 
on particle configuration at all and results in the equiprobable ensemble. 
Apparently the ZRP under consideration is not the case like this. The 
way how to reduce the problem to the one with uniform groundstate was 
proposed in The main idea is to look for the solution in the form, 

Ft(C7)=P,t(C)i^°(C7). (24) 

Here Pst (C) is the stationary probability defined in Q. It is not difficult 
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to check that Fj^ (C) satisfy the foUowing equation 

N 

xF° {ni - fci + ^2, . . . ,nAr - fcjv + fci) . (25) 

This form of the equation has two important advantages. The first is 
that the coefficient before F^ (C) under the sum is the product of one 
site factors. The second is that if we formally define p{Q) = 0, the values 
of ki-s will take on values of and 1 with no respect to the value of rii. 
To proceed further we should go to a different way of representation of 
system configurations. Let us define the configuration by the coordinates 
of particles written in non-decreasing order 

C = {xi, . . . , xi < X2 < ■ ■ ■ < XM- (26) 

Obviously this is nothing but the change of notations, both representations 
being completely equivalent. Below we refer to them as to the occupation 
number representation and the coordinate representation respectively. We 
are going to show that the eigenfunctions of the equation (|25|) written in 
the coordinate representation can be found in form of the Bethe ansatz. 
We first consider the cases M = 1 and M = 2 subsequently generalizing 
them for an arbitrary number of particles. 

A. The case M = 1. 

The eigenfunction problem for the master equation for one particle is 
simply the one for the asymmetric random walk in discrete time. 

A (7) F° (x) = e>F° (a; - 1) + (1 - p) F° (x) (27) 

The eigenfunction is to be looked for in the form 

F°(x) = z-^ (28) 

where z is some complex number. Substituting H28|l to l|27fl we obtain the 
expression of the eigenvalue 

A (7) = e^pz + {l-p). 

The periodic boundary conditions 

F° (x + N) = F° (x) 

imply the limitations on the parameter z 
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B. The case M = 2. 

In this case we consider two particles with coordinates Xi and X2- If 
Xi ^ X2 the equation is that for the asymmetric random walk for two 
non-interacting particles 

A(7)F°(xi,X2) = (e»'F"(xi-l,a;2-l) 

+e> (1 - p) {F° {xi -l,X2)+F° {xi,X2 - 1)) 

+ {l-pfF'^ixi,X2), (29) 

while the case xi = X2 = x should be treated separately. 

A (7) F° {x, x) = e> (2) F° (x - 1, x) + (1 - j5 (2)) F" {x, x) (30) 

The general strategy of the Bethe ansatz solution is as follows. We want to 
limit ourselves by the only functional form of the equation for F° (C) of the 
form (|29|) . To get rid of the additional constraints imposed by interaction, 
like that in H30I) . we recall that the function F° {xi, . . . ,xm), is defined in 
the domain xi < X2 < • ■ ■ < xm- If we formally set xi = X2 = x in the 
equation (|29|l . one of the terms in r.h.s., F° {x,x — 1), will be outside of 
this domain. We could redefine it in such a way, that the Eq. (|Sn|) would be 
satisfied. This leads us to the following constraint. 

aF° {x, x) + hF° (x, x - I) + cF^ {x -l,x) + dF° (a; - 1, a; - 1) 0, (31) 

where 

« = ((1-P)'-(1-P(2))), 

b = eXl -p) , 

c = e^(p(l-p)-p(2)), 

d = {e^pf 

The solution of the equation (|29|) can be looked for in the form of the Bethe 
ansatz, 

F° (XI,X2) = Ai,2Zr'^2""' +^2.l2r'22""'- (32) 

Substituting it into Ea. (|29|l we obtain the expression for the eigenvalue 

A (7) = (e>zi + (1 - p)) (e>Z2 + (1 - p)) , 

while the constraint (|31|l results in the relation between the amplitudes 
Ai^2 and ^2,1. 

^1,2 _ a + bzi + CZ2 + dziZ2 
^2,1 a + bz2 + czi + ^2:12:2 
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The last step, which makes the scheme self-consistent, is to impose periodic 
boundary conditions 

which lead us to the system of two algebraic equations. The first one is 

_ a + bzi + CZ2 + dziZ2 
^ a + bz2 + czi + dziZ2 

and the second is obtained by the change zi < — > Z2. 

C. The case of arbitrary M. 

For an arbitrary number of particles we will follow the same strategy. 
We consider the equation for noninteracting particles, which jump forward 
with the probability p 

M 

a{^)f°{xi,...,xm) = En(e»''(i-p)'"'' 

xFO(a;i-fci,...,XM-fcAf). (33) 

Here, all the numbers fc^-s run over values 1 and 0. Our aim is to redefine 
the terms beyond the physical domain, (|26l) . in terms of ones within this 
domain in order to satisfy the equation with ZRP interaction (|25|) written in 
the coordinate representation. This redefinition results in many constraints 
on the terms (•■•)■ it will be seen below, the Bethe ansatz to be 
applicable all these constraints should be reducible to the only one, 

= aF°{...,x,x,...) + bF°{...,x,x-l,...) 

+cF" (. . . , a; - 1, x . . .) + dF" (. . . , a; - 1, X - 1, . . .) , (34) 

which has been found for the two particle case. To this end, we first recall 
that the r.h.s. of Ea. (|25|l is the sum of the terms F^ (. . .) corresponding 
to configurations, from which the system can come to the configuration 
in l.h.s. with coefficients factorized into the product of one-site terms. 
Therefore, the processes corresponding to a particle arriving at a site can 
be considered for each site separately. Consider, for instance, Ea. H25|) . 
where the argument of the function in l.h.s. is a configuration with a site 
X occupied by n particles. 

A (7) (. . . , (xf ,...)- E{,, J ni^. ("^))'' (^-p 

X [e>(n)F«(...,(x-l),(x)"-\...) (35) 
+ (l-p(n))F°(...,(x)",...)] 
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The terms of r.h.s. can be grouped in pairs shown in square brackets, which 
correspond to the processes, in which a particle comes or does not come 
to the site x. Here, (x)" denotes n successive arguments equal to x, i.e. 
the site x is occupied by n particles. The primed summation and product 
run over all sites apart of x. Obviously the coefficients in the equations 
for noninteracting particles can be factorized in similar way, such that the 
term corresponding to the site x looks as follows 



ki=0 fc„=o 



(36) 



x^0(. 



ki , . . . , X ,...,). 



We want this term to be equal to the one in the square brackets of Ea. H35|l . 
If n = 2 the form of the term in square brackets is similar to the r.h.s 
of Ea. (|30|l and can be treated (i.e. reduced to the noninteracting form 
(|36|l ) using the constraint 1)34(1 . The equality for general n can be proved 
by induction. Suppose it is valid for n — 1. Then we can perform the 
summations in Ea. (|36|l over fc^-s for i ~ 2, . . . resulting in 



(1 - p) e> {n - 1) F° . . , (x - 1) , (x) 
+ (l-p(n-l))F°(...,(x)",...)] 



n-2 



(37) 



+e> [e> (n - 1) F" (. . . , (x ~ 1)' , (x)""' , . . .) 
+ (l-p(n-l))FO (...,(x-l),(x)"-\...)" . 

The first summand contains the term F'^ ^. . . , x, (x — 1) , (x)" ^ , . 

which, being beyond the physical region (|26() , can be expressed in terms of 
ones inside the physical region by using the constraint H34|l . Equating the 
coefficients of the terms F'^ (. . .) with the same arguments of Eas. (|36|) and 
H37|l we obtain the relation between p (n) and p{n ~ 1) 



p{n)=p + qp{n-l), 
where we introduce the notation q defined as follows 

p{2)=px{l + q). 



(38) 



(39) 



The recurrent relation (|38|l can be solved in terms of q— numbers as was 
claimed in Eas. (|9ll0|l . 

Thus, we have shown that in the physical domain (|26|l the solution of 
the equation for ZRP coincides with the solution of the equation for non- 
interacting particles if the latter satisfies the constraint (|34|l . Then we can 
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use the Bethe ansatz 

M 

F^{x,,...,xm)= J2 ^-i.-.-Mll^-r (40) 

{o-i,...,ctm} *=1 

for the eigenfunction of the free equation. Here zi, . . . ,zm are complex 
numbers, the summation is taken over all p\ permutations {cti, . . . jCr^/} 
of the numbers (1,...,M). Substituting the Bethe ansatz H4U|I into the 
Eq. H33|) , we obtain the expression for the eigenvalue 

M 

Ah)^l[{e^pz^ + {l-p)). (41) 

i=l 

Substituting it to the constraint H34() we obtain the relation between pairs 
of the amplitudes ACTi,...,CTj,f , which differ from each other only in two indices 
permuted 

A...ji... a + bzj + czi + dZiZj 

— = Si. = ; — -.. (42) 

A...i^j... a + bzi + czj + aziZj 

With the help of this relation one can obtain all the amplitudes ^cri,...,o-M 
in terms of only one, say Ai^...^m, by successive permutations of indices, 
which results in multiplication by the factors Sij. For example, for three 
particle case we have 

^2,1,3 = "Si, 2^1, 2, 3, ^2,3,1 — '5'l,3'S'l,2^1,2,3, ^3,1,2 — '5'2,3'S'l,3^1,2,3 
^1,3,2 = 'S'2,3^1,2,31 ^3,2,1 = '5'2,35'i, 35*1,2^1, 2, 3- (43) 

In general, the procedure to be uniquely defined, it should be consistent 
with the structure of the permutation group. Specifically, if Bi is an el- 
ementary transposition that permutes the numbers at i-th and {i + l)-th 
positions, it satisfies the following relations. 

^\ - 1 (44) 

If the numbers at positions i, (i + 1), and (i + 2) are j, /c, and I respectively, 
two relations of Eq. (|44ll are reduced to 

Sj.kSj.iSkd = Sk,iSj.iSj^k, 

Sj,kSk,j = 1, (45) 

which are apparently true. 

The last step we need to do is to impose the periodic boundary 
conditions. 

F\xi„...,xm) =F°{x2,,...,XM,xi+N) (46) 
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They are equivalent to the foUowing relations for the amplitudes. 

A -A 

-^ai,...,crM ~ ^o'2,...,(Tm,o-i^cti 

Consistency of this relation with Eq. H34|) results in the system of M alge- 
braic Bethe equations, 

^-jv _ ^ ^M-iYi a + bzi + cZj+dziZj ^^^^ 
' a + bz^ + czi + dziZj 

We should emphasize that the crucial step for the Bethe ansatz solvability 
is the proof that all many-particle interactions can be reduced to the two- 
particle constraint, Ea. H34|) . Existence of any other constraints on the 
eigenfunctions would result in new equations for parameters Zi , which would 
make the resulting system of algebraic equations overdetermined. 



D. Bethe ansatz for the ASEP. 



As it was noted above, the Bethe ansatz solution for the ASEP is quite 
similar to that for the ZRP and can be done in parallel. Indeed, if we 
define the function iC) in the same way as it was defined for ZRP, the 
equations for it can be obtained from those for ZRP by the variable change 



{xiX2, ■ ■ - .Xm} {XI,X2 + l,...,XM+ M - 1} , (48) 

which corresponds to the ZRP-ASEP transformation described above. 
Thus, the free equation does not change its form, with the only difference 
that the physical domain for the ASEP implies every site to be occupied 
at most by one particle. 

Xi < X2 < ■ ■ ■ < XM- (49) 

Therefore the Bethe ansatz (|40|l substituted to the free equation results in 
the same form of the eigenvalue H41|) . The two particle constraint, which is 
used now to redefine the nonphysical terms containing the pair {x,x), have 
the following form 

= aF° {x, x + l) + bF° {x, x) + cF° {x-l,x+l} + dF" (x -l,x). (50) 

If we insert here the Bethe ansatz (|4Uf) . we obtain the following relation for 
the amplitudes 
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All the other arguments, which extend the problem to the general 

M— particle, case arc completely the same as those for the ZRP. Finally, 
the periodic boundary conditions lead to the Bethe equations 

^ az^ +b + cziZj + dzi 

Here we recall that the lattice size is equal now to L = A'' + M rather than 
N. This equations can be rewritten in the following form 



a + hzj + czi + dziZ 

where we introduce the notation 

M 



T=n 



(54) 



This term is the only difference between the Bethe equations for the ZRP 
and the ASEP. Taking a product of all the equations we obtain, 

= 1. (55) 

The term T has a very simple physical meaning. This is the factor that cor- 
responding eigenfunction multiplies by under the unit translation. Thus, 
the eigenfunction that corresponds to the set {zj} satisfying a relation 
T = 1 is invariant with respect to any translations for both ZRP and ASEP. 
This is what we meant mentioning the translational invariance above. Ap- 
parently the solutions of the Bethe equations, which satisfy the relation 
T = 1 coincide for the ZRP and the ASEP and so do any quantities con- 
structed with them. Below we consider such quantity and, therefore, do 
not make a distinction between the ZRP and the ASEP. 



IV. THE LARGEST EIGENVALUE. 

To proceed further, we introduce new variables yi defined as follows 

(56) 
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where A is defined in Ea. H14|) . In these variables the Bethe equations (|47|l 
and the eigenvalue (|41() simplify to the following form 

, \ -AT 1\I 



i + xyj ^ ^ ^ l\y,-m' ^^^^ 

A (7) = (58) 



i=l 



Let us now consider the eigenstate, which corresponds to the maximal 
eigenvalue Aq (7). Note that in the limit 7 — > the equation (|7J) for Ft{C) 
turns into the Markov equation for the probability |0J. The largest eigen- 
value of the Markov equation is equal to 1. The corresponding eigenstate 
is the stationary state |SJl- It then follows from the definition, that F^{C) 
in this limit becomes constant, i.e. it is the same for all configurations, 
and the corresponding solution of the Bethe equations is Zi — I or yi = 0. 
If 7 deviates from the analyticity and no-crossing of the eigenvalue is 
guaranteed by Perron-Frobenius theorem. By continuity we also conclude 
that the parameter T defined above is equal to 1 for arbitrary 7, 

. J; 1 -I- Ay, 



A. The case q — 0. 

In the case q — 0, the form of Eg. (15 711 allows one to obtain the largest 
eigenvalue for arbitrary M and TV [13. If we define the parameter, 

M 

B = i-f-'e^^l[y,, (60) 

the solution of the Bethe equations, which corresponds to the largest eigen- 
value, will be given by M roots of the polynomial 

P{y) = {l + Xyf B~{l-yfy'\ (61) 

which approach zero when B tends to zero. Then the sum of a function 
/ (x) analytic in some vicinity of zero over the Bethe roots, can be calcu- 
lated as the following integral over the contour closed around the roots: 
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Particularly, after the integration by parts the expression for In Aq (7) has 
the following form: 

lnAo(,)^/i^-Vlnfl-^(^+y"V (63) 

At the same time the expression for 7 as a function of B can be obtained 
by taking a logarithm of Eq. H59|l 

^ I dy f 1 A \ / B(l + Xy)'^\ 

■'-Mf27,[—^—yy'['-jr-ir^)- 

To evaluate it we make a series expansion of the logarithms in powers of 
B and integrate the resulting series term by term. The procedure is valid 
until the resulting series are convergent. The resulting series, obtained after 
some algebra by using standard identities for hypergeometric functions, 
are given in Eas. H15ll6|) . The large time asymptotics of the cumulants 
of the distance travelled by particles limt^oo {Y")^/t can be obtained by 
eliminating the parameter B between two series. For instance, the exact 
value of the average flow of particles is of the following form: 

(1-M.l-N,_ 

N t-oo t L-1 ( -M, -N ^ ' 

1^ 1-L 

In the thermodynamic limit N — > 00, Af — > oo,M/L ~ c the average 
velocity of particles v = N(t)/M saturates to the formula obtained in 0, 



, = ^^^l^^IE^^EA. (66) 

M 2c ^ ' 

It is straightforward to obtain the other cumulants up to any arbitrary 
order. It is not difficult also to study the asymptotic behavior of the series 
in the large N limit. To this end one can evaluate the saddle point asymp- 
totics of the integrals instead of their exact values or use the asymptotic 
formulas for the binomial coefficients and the hypergeometric functions. 
The thermodynamic limit however allows significant simplification already 
at the stage of writing of the Bethe equations, which makes possible to 
consider more general case of arbitrary value of q. 



B. The case of arbitrary q. 



The technique used in this section was first developed to study the non- 
factorizable Bethe equations for the partially ASEP [13 . IT^ and then ap- 
plied to the study of the ASAP [s^- Since the analysis is quite similar. 
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we outline only the main points of the solution. Technical details can be 
found in the original papers. 

The scheme consists of the following steps. We assume that in the ther- 
modynamic limit, N oo, M oo, M/N — p, the roots of the Bethe 
equations (|57|l are distributed in the complex plain along some continuous 
contour F with the analytical density R{x), such that the sum of values of 
an analytic function f(a;) at the roots is given by 

P r 
Vf(a;,) = N f{x)R{x)dx. (67) 
z=i -'r 

After taking the logarithm and replacing the sum by the integral, the sys- 
tem of equations (|57|l can be reformulated in terms of a single integral 
equation for the density. 

- In (i^) + 7 - / In (^) Riy)dy = - 2.^ £ R (x) (68) 

The r.h.s. comes from the choice of branches of the logarithm, which 
specifies particular solution corresponding to the largest eigenvalue. The 
choice coincides with that appeared in the solution of the asymmetric six- 
vertex model at the conical point and later of the ASEP It can 
be also justified by considering the limit q — 0, where the locus of the roots 
can be found explicitly. The point xq corresponding to the starting point 
of roots counting can be any point of the contour. The integral equation 
should be solved for a particular form of the contour, which is not known 
a priory, and being first assumed should be self-consistently checked after 
the solution has been obtained. In practice, however, an analytic solution 
is possible in the very limited number of cases. Fortunately, we can proceed 
by analogy with the solution of the asymmetric six-vertex model at the 
conical point 33] choosing the contour closed around zero. The solution of 
the equation H68|l yields the density 

R^"\x)^^{p-g,,,{x)), (69) 
where the function gq \{x) is defined as follows 



This case corresponds to 7 = and hence Ao(7) — 1. Since R^'^^x) is 
analytic in the ring < |a;| < A, the integration of it along any contour 
closed in this ring does not depend on its form. Therefore, to fix the form 
of the contour additional constraints are necessary. Such a constraint was 
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found by Bukman and Shore |33 | as follows. They assumed that when 7 
deviates from zero, the contour becomes discontinuous at some point Xc- 
It is possible then to solve Eq. (|68|l perturbatively considering the length of 
the gap of the contour as a small parameter. It turns out that the solution 
exists only if the break point Xc satisfies the equation 

R'^°\x,)^Q, (71) 

which fixes the location of Xc as well as the location of all the contour T in 
the limit 7^0. This method however allows a calculation of only the first 
derivative of Aq (7) at 7 = in the thermodynamic limit, i.e. the leading 
asymptotics of the average flow of particles. To study the behaviour of 
Ao(7) for nonzero values of 7, a calculation of the finite size corrections 
to the above expression of {x) is necessary. To this end, we use the 
method developed in [T3| . In fact, once we have the expression for 
R^^\x), Ea. (|5^ . it can be directly substituted to the formula for the finite 
size asymptotic expansion of R{x) obtained in 32] for the ASAP. As a 
result we obtain the following parametric dependence of R (x) on 7, both 
being represented as the functions of the same parameter k, 



Rs = 



X 

n=0 



_l 1 

iV3/2 2^1-9l«l 
/ r(n+ I) C2n+l,^ 



y[ — ] -V"-"^2y^2n+l,. . ) ( ^ 

_L ^/ i V r(n+|)c2«+l T- / N /VON 

7 = -^21^[2n) ,n+f V2l ^^"+^^"^ ^^^^ 
n— ^ 

Here Lia(x) is the polylogarithm function, 

00 g 

Lia(x) 



s—1 



Rs and Ri^^ are the Laurent coefficients of R{x) and R^-^'^{x) respectively 
defined as follows 

R{x)= Rs/x^+\ (74) 



s— — 00 



and Cn,s and c„ are the coefficients of in {^^^^akX^Y and 
log (^^Q afex'^) respectively, where a„ are the coefficients of the inverse 
expansion Z~^{x) of the function Z{x) ^ f^^R{x)dx near the point 
Z{xc). As the derivative of Z{x) vanishes at a; = Xc, Ea. H71|l . the inverse 
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expansion is in powers of \J x — Z{xc), 

oo 

Z~^{x) = ^ a„ (a; - Z(a;c))"^^ . 

Being substituted to the Abel-Plana formula 12], which is used to evaluate 
the difference between the integral in Ea. H68(l and the sum in the original 
Bethe equations, it becomes a source of 1/VN corrections to The 
location of Xc is to be self-consistently defined from the equation R (xc) = 
0. For the first three orders of expansion the coefficients a„ can 

be obtained from the inverse expansion of zero order function Z(°)(a;) = 
— /^^ {x)dx, while 7?^°^ (x) has been obtained above. The details of 
calculations can be found in [s^- The sum over the roots of a function 
f(2;) analytic inside the contour can be evaluated in terms of the Laurent 
coefficients Rs of R (x) 

M oo 



^f(a;,) = 27:iNj2l''isRs- (75) 

i=l s=l 

and the of Taylor coefficients of i{x) 

oo 

f(.T) = ^f,x^ (76) 

s=l 

This allows evaluation of In Aq (7). Finally the point Xc enters into all the 
results through the coefficients C2n+i.s and C2n+i- The final expression for 
Aq (7) obtained in the scaling limit jN^^^ — const has the form (|17I19() . 
which was obtained earlier for the ASEP and the ASAP and claimed to 
be universal for the KPZ universality class The constants 0, fci, k2 

specific for the model under consideration. 



\Xr 



1 + XXc 



(77) 



Xc 2^ 9',,x(.Xc) + {l + Xxc)Xg'lM 



ki = -^^^^ ^ --^-^ (78) 

(l + AXe)^ff;,(Xe))5/2 

^2 = ^2TTXcgg x i^c), (79) 

are expressed through the derivatives of the function gq x{x), fEa. (|70|l ). As 
it follows from the explicit form of Ea. (|71|l . the unphysical parameter Xc is 
related to the density p by the equation. 



(80) 
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In Sec. VI we show that relations (|77I80|1 can be obtained from the par- 
tition function formahsm for the stationary state without going into the 
Bethe ansatz solution. The parameter Xc then turns out to be related to 
the fugacity of a particle in the infinite system. Thus, this result follows 
from the properties of the stationary state only, rather than from the solu- 
tion of the dynamical problem. In general the partition function formalism 
allows the calculation of any stationary spatial correlations. The Bethe 
ansatz, however, allows one to probe into really dynamical characteristics, 
i.e. those, which are related to unequal time correlation functions, such 
as higher cumulants of the distance travelled by particles. They can be 
expressed through the constants fci, fc2 from Eas. (|78l79|l by taking deriva- 
tives of the formula (|15|l . which yields the following behavior, specific for 
the KPZ universality class, 



Studying the dependence of the results on the physical parameters p and 
q one should solve the equation H8U|) to find the behavior of Xc- Even 
without the explicit solution we can say that for generic values of p and 
q the point x^ is located at the positive part of the real axis in the seg- 
ment [0,1), where the function gq \ {x) increases monotonously from to 
infinity. Apparently the constants 0, ki , /c2 do not have any singularities 
for Xc varying in this region and, thus, no abrupt change of the behavior is 
expected. Particularly, the results obtained reproduce the corresponding 
results for several models studied before. For example, in the continuous 
time limit, p — > 0, keeping only the first order in p we obtain the corre- 
sponding constants for g— boson asymmetric diffusion model |l7j|.|27j. which 
itself contains the drop push model and phase model |28l| as particular 
cases, and can be related also to the totally ASEP. The only interesting 
exception is the limiting case p — > 1, which turns out to exhibit a kind of 
phase transition phenomena. 

V. ASYMMETRIC AVALANCHE PROCESS AND PHASE 
SEPARATION IN THE DETERMINISTIC LIMIT. 

Let us first consider qualitatively what happens in the limit 



p = 1 — St, St — > 0, 



(81) 



in the domain 



1 < g < 0, p < 1. 



(82) 



In this limit a particle from a single particle (SP) site, i.e. from a site 
occupied by one particle only, almost definitely takes a step forward at 
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every time step. At the same time, one particle from a many particle 
(MP) site, i.e. from a site occupied by more than one particle, takes a 
step with the probability in general less than one. When the jumps of 
particles from SP sites are purely deterministic, i.e. the equality p = 1 
holds exactly, the ZRP dynamics becomes "frozen" as soon as the system 
arrives at any configuration consisting of SP sites only. By "frozen" we 
mean that at every time step particles from all sites synchronously jump one 
step forward. Thus, the structure of the configuration remains unchanged 
up to a uniform translation. Therefore such configurations play the role of 
absorbing states when p — I. When p is less than 1 by a small value 6t, 
the system can go from one absorbing state to another with the probability 
of order of St. This happens if at least one particle decides not to jump. 

It is easy to see that the limit under consideration is directly related to 
the ASAP. Let us look at the system from the comoving reference frame, 
which moves a step forward at every time step. In this frame particles 
from the SP sites either stay with the probability (1 ~ 5t) or take a step 
backward with the probability 6t. At the same time the MP sites play a 
role of columns (or avalanches) of particles moving backward. Individual 
history of every such a column develops according to the following rule. 
If n > 1 particles occupy a site, either all n particles jump to the next 
site backward with probability = (1 — p (n)) or {n — 1) particles jump 
and one stays with probability (1 — /i„). One can see that the definition of 
probabilities /i„ coincides with those for the ASAP, Ea. (|13|l . in the leading 
order in 6t. If we associate one step of the discrete time t with the con- 
tinuous time interval 5t, neglect the terms smaller than 6t in the master 
equation and take the limit 5t ^ 0, we come to the continuous time master 
equation for the ASAP . Note, that in the definition of the ASAP at 
most one avalanche is allowed to exist at a time. In principle, the discrete 
time dynamics allows for more than one simultaneous MP sites. However, 
the probability of formation of two or more of them for a finite number of 
steps is of order of (St)'^ and, as such, corresponding processes vanish when 
the limit ^ is taken. 

As shown in the ASAP exhibits the transition from the intermittent 
to continuous flow regime. In this transition the average avalanche size (the 
number of particle jumps in an avalanche) diverges in the thermodynamic 
limit when the density of particles approaches the critical point from below. 
To explain what this behaviour means in terms of the discrete time model 
under consideration, we note that the average avalanche size is roughly 
equal to the average number of particles occupying an MP site times the 
average number of steps, which an MP site persists for. Let us consider 
a configuration with one MP site on the lattice. In average, the flow of 
particles into the MP site is equal to the density of mobile (solitary) parti- 
cles. The flow out of the MP site with n particles is p{n), which saturates 
exponentially rapidly to the limiting value p{oo) = 1/(1 — g) when n grows. 
Therefore, if p < p^, = 1/(1 — g), the outflow exceeds the inflow and the 
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FIG. 2: The flow-density plots for q = -1/2 and 5t = 0.5,0.5 x 10"\0.5 x 
10"^0.5 X 10-^0.5 X 10"'', 0.5 x 10"^0, going in bottom-up direction respec- 
tively. 

MP site tends to decay for a few steps. As a result the average occupation 
number and the life time of an MP site below the critical density are finite. 
On the contrary, if p > p^, the outflow is less than inflow and the dynamics 
favours the growth of the MP sites. On the finite lattice the MP site grows 
at the expense of the density of solitary particles until it absorbs enough 
particles to equalize inflow and outflow. This happens when the MP site 
contains N{p — p^) particles. Any deviation from this number destroys the 
balance between the inflow and outflow and the system tends to restore 
this value again. As a consequence the life time of such an MP site grows 
exponentially with the system size "351. 




Of course the existence of only one MP site on the lattice is specific for 
the continuous time ASAP, where the probability of appearance of two MP 
sites simultaneously is neglible. In the discrete time case a possibility to 
have many MP sites at the lattice should be considered as well. However, 
if the characteristic occupation number of MP sites is such a large that the 
outflow from MP sites can be roughly replaced by p{oo), the above claim 
based on the balance between inflow and outflow remains valid: the total 
number of particles in all MP sites is equal to N{p — p^). As a result, 
the particle current growing linearly with the density below the critical 
point, stops growing at p^ due to appearance of MP sites, the occupation 
numbers of which grow instead. Of course these arguments describe the 
limiting situation p — s- 1. When St is finite, the transition is smoothened. 
As shown in Fig.(|21), the flow-density plot, obtained by substituting the 
result of numerical solution of Ea. H8U|) to Eq.||77)), approaches the form of 
broken line consisting of two straight segments: (p = p for p < p^ and 
= Pc for P > Pc- 

Of course the above qualitative arguments being far from rigorous can 
only serve as an illustration. A specific distribution of the occupation num- 
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bers of MP sites can be obtained with the help of the canonical partition 
function formalism, which is described in details in the next section. In 
present section we analyse the behaviour of generating function Ao(7) of 
the cumulants of the integrated particle current Yt that can be extracted 
form the solution of the Bethe ansatz equations discussed above. 

To this end, we note that the universal scaling form (|17I19|I of Ao(7) 
holds in all the phase space, provided that the parameters A:i,fc2 are 
finite. What depends on the parameters of the problem is the typical scales 
encoded in (j), ki,k2, which characterize the average and the fluctuations of 
the integrated particle current Yt. 

Before going to the analysis let us look how the results are related to 
those for the ASAP. To this end, we recall that in order to transform the 
ZRP to the ASAP we consider comoving reference frame, which takes one 
step forward at each time step. The distances, which particles travel in 
these reference frames, are related to each other as follows 

r/^^ = Aft-r,^^^^. (83) 

Hence, the largest eigenvalues of the equations for the generating functions 
Ft (C) for the ASAP and ZRP are related as follows 

lnAo(7) =7M-'5rA^^^^ + 0(<5T2). (84) 

By A^'^'^'^ we mean the largest eigenvalue of the continuous time equation 
for the ASAP obtained in 32] . The factor St which comes with A^'^'^'^ 
reflects the time rescaling, in which one discrete time step is identified with 
the infinitesimal interval St. The formula (|84ll suggests that the expansion 
of In Aq (7) in powers of St exists, which turns out to be true only in the 
subcritical region p < p^. In what follows we solve the equation (|5n|l in 
the limit p — > 1, i.e. 6t — > 0. Once the value of Xc is obtained from this 
solution, it is to be directly substituted to the formulas H77I79|I to obtain 
</>, ki,k2. 

When St ^ the function gq \ (x) can be summed up to the following 
simple form 

X 

9q-l/q{x) = , (85) 

X — q 



such that Ea. (|80|l can be easily solved 

qp 



p-1 



(86) 



Then Eas. H77l79ll yield: (j) = p and ki = 0, i.e. IuAq = jM. This corre- 
sponds to synchronous jumps of all particles without any stochasticity. For 
small nonzero St we expect the correction of order of St to appear. This 
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should be true at least if the argument of gq^\ (x) is away from any of its 
singularities, which turns out to be the case only in the subcritical region 
p < p^. In general in this case one can represent x (x) as a Taylor series 
in powers of 5t with non-singular series coefficients. Then Ea. H8U|l can be 
solved perturbatively order by order in 5t. To this end, we look for the 
solution in the form of the perturbative expansion in powers of 5t 



The first order solution of Eq. H8U|) yields 



(87) 



(1-9)P 



1 



which results in the following values of 0, ki,k2'- 



4) ~ p + St 



St- 



(1-g) 

(9-1) 



IP 



- i-g^ \p-ij ^ 

s^{s -l + 2p) ( qp 



8^<z((l-p)p) 



3/2 



k2 ~ y/2TT{l - p)p. 



(89) 



(90) 
(91) 



Thus, the first order calculation allows taking into account the fluctuations 
of the particle current due to spontaneous MP site formation. The values 
obtained indeed reproduce from Js^l . Close to the critical point x^c^ 

diverges, 

X^c ^ ^ {Pc - PV"^ > 

resulting in the divergent contributions to 0, fci, /c2 

^^''^ - {Pc - pv^ , fci'^ - {Pc - pv' , fcf^ - iPc - py' ■ 

The reason of these divergencies have a clear physical meaning. For ex- 
ample, the first order correction to the flow (j) gives a fraction of particles, 
which stay in MP sites. Roughly speaking it is a product, tMpnMPPMP, of 
the life time of an MP site Imp, its average occupation number ump and 
the probability of appearance pmp, the latter being of order of 5t. As it 
follows from the above arguments based on the balance between the flows 
into and out of an MP site, the two former. Imp and ump, should diverge 
as the density of particles approaches its critical value. 

The singularities in the solution limit the applicability of the perturbative 
scheme used. It to be applicable the solution Xc should be far enough from 
the first singularity oi gq^\ (x), at the positive half-axis, a; = 1. Specifically, 
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the parameter (^tx^c^ / {p^ — p)j should be small to ensure the existence 
of nonsingular expansion for gq x (x). This yields 

St{p^ - « 1, (92) 

resulting in 

l-Xe~/9, -p>^t1/^ (93) 

To proceed further, we note that the function gq^\ {x) monotonously 
grows along the positive half-axis from zero at the origin, a; = 0, to infinity 
at the pole x = 1. As the function gq^\ (x) runs over all positive real values 
between these two points, including the value of p, the solution of Ea. l|5n|l . 
Xc, is always located at the segment < x < 1. According to the above 
perturbative solution almost all this segment is exhausted by the subcritical 
domain, Ea. (|93l) . except the small vicinity of x = 1, i.e. (1 — Xc) ^ St^^^. 
Therefore, this vicinity is where Xc should be looked for at the critical point 
and above. 

This remark can be used as a basis for another perturbative scheme. Let 
us suppose that 

= 1 - A, (94) 

where A is a small parameter. The function gq^\ (x) can be represented in 
the following form 

fc=0 ^ k=0 ^ 

This is obtained by expanding the expression under the sum in Ea. l|7U|) in 
powers of g" and then changing the order of summations. In this form 
gq^x (x) has a form of the sum of terms, each with one of the following 
simple poles 

x',=q-^ 4! = -<7-*A-\ fc = 0,1,2,... 

The pairs of the poles x'j. and x'l_-^ merge when 6t tends to zero, so that 
corresponding terms having equal absolute values and opposite signs cancel 
each other. As a result, for finite A the only term of gq \ (xc), which 
survives in the limit 6t 0, is the one with the pole Xq (see Ea. (l85|) '). On 
the other hand the term under the first sum in Ea. H95|l with the pole Xq = 1 
tends to infinity when A tends to zero. The compromise is to consider the 
two limits simultaneously. Then the divergent contribution from the pole 
x'q can be cancelled up to a finite constant by the term with the pole 
x'l — — (qX) . The value of the resulting constant, which depends on the 
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relation between St and A, should be chosen such as to satisfy Ea. H8U|) . 
Substitution of Ea. (|94|l to Ea. HSOl) shows that the value of gq^\ (xc) to be 
finite in the limit A ^ 0, the following ratio should be kept constant. 

St/A'^ const (96) 

Using this fact we can look for the relation between St and A in form of 
the expansion in A 

St^{1- l/q)-^ (aA2 + /3A3 + 0(A*)) . (97) 

The coefficient (1 — is for the further convenience. Solving the equa- 

tion (|80|l order by order in powers of A we fix a and /3, 

a — p — Pc (98) 
(3 = a{l + a)-q{l-qy\ (99) 

Inverting the relation H97|) we express Xc in terms of St. Finally we substi- 
tute Xc into Eas. H77l79() . The cases a = and a > should be considered 
separately. 

- In the case a > 0, which corresponds to the density above critical, 
p > p^, we obtain in the leading order 



Pc 



St^I"" 



2n \ni/2 



Pc 



(100) 



fci ^^Jl-pj3/4^rV4 ^ ^^^^^ 

fc2 ^ 2^^F(l-pJl/^i^^-|^. (102) 

As follows form these formulas the flow of particles 4> in the thermodynamic 
limit is equal to the critical density p^ as expected. It is possible also to 
obtain correction to the thermodynamic value of the particle flow 
which was shown to be universal being proportional to the nonlinear 
coefficient A of the KPZ equation |3^ . 

1 3 pjl - pj 

The value of the correction diverges in the critical point. Surprisingly 
its value does not depend on St. Higher cumulants, responsible for the 
fluctuations around the average flow grow when St tends to zero, 

ivnx (n^ n ^(3n-7)/4 
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- In the critical point a — Q, i.e. p — p^, we obtain, 

^ P,-5t'''[pI{1-P,)]"\ (105) 

kr (|£)'^'[547^(l-pJpJ-l/^ (106) 

fc2 ^ V6^(l-pjp,. (107) 
The 1/A^ correction to the average flow looks as follows 



1 1 /p,xV' 



(108) 



growing with the decrease of 5t. The other cumulants have similar 5t 
behavior, 

t— ►oo t 

As follows from the results obtained, though the average flow saturates to 
a constant value at the critical point, its fluctuations grow when 5t goes 
to zero. 

The generating function of cumulants An (7) can be related to the large 
deviation function for the particle current [l3|. As a result the probability 
distribution of Yf /t reads as foUows 



P ( y = y ) = exp 



tc \ G 



(110) 



where the universal function H{x) is given by the following parametric 
expression 

H{x) = pG'{P)~G{(3) (111) 
X = G'{I3). (112) 

with two model dependent constants 

t^ = A^3/2^-l 

G = kik2. 

The function H{x) behaves as (a; — 1)^ for x ^ 1 and as x^^^ and x^/^ as 
X goes to plus and minus infinity respectively. One can see from Ea. HllU|l 
that the parameter G plays the role of the characteristic scale at which 
the current fluctuations become non-Gaussian. This takes place when the 
argument of H{x) in Ea. H110|l becomes of order of 1. Note that G coincides 
with 1/A^ correction to the average flow. This is a reflection of the fact 
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that the correlations between particle jumps, which cause deviations from 
the Gaussian behaviour, owe to the finiteness of the system and to the 
periodic boundary conditions. The parameter tc is the characteristic time 
in which such fluctuations become unlikely. The characteristic time tc being 
of order of iV^/^ signifies the KPZ behaviour characterized by the dynamical 
exponent z = 3/2. In fact, for the totally ^ and the partially JJj ASEP 
the inverse time coincides with the next to the largest eigenvalue of the 
master equation up to a constant of order of 1. We expect this to hold also 
in our case. Below the critical point p < we have tc ~ N^^'^6t~^ and 
G ^ St. Above the critical point, p > Pc, the characteristic time scales as 
tc ^ N'^/'^St^'^/^ growing with the decay of 6t and the fluctuation scale is 
finite G ^ {p — Pc) ^ growing as the density approaches the critical point. 
Exactly at the critical point we obtain G ^ 5t-^'^ and tc ~ St^'^N^I'^. 
One can see that below the critical point even very small fluctuations are 
non-Gaussian. Above the critical point the range of Gaussian fluctuations 
is finite. At the critical point the fluctuations remain Gaussian in a very 
wide range. On the other hand the time of decay of the fluctuations at 
the critical point is much smaller than below and above. 



VI. CANONICAL PARTITION FUNCTION FORMALISM FOR 
THE STATIONARY STATE. 

In this section we show that, the partition function formalism |7|.[38| 
allows a calculation of some physical quantities yielding the same results 
as obtained above. Particularly the results obtained in the saddle point ap- 
proximation are equivalent to the results obtained from the above solution 
of the Bethe equations in the thermodynamic limit. We also analyze the 
change of the occupation number distribution under the phase transition 
in the limit 5t ^ {). 



A. Canonical partition function and stationary correlations. 

The partition function of the ZRP is defined as the normalization con- 
stant of the stationary distribution ||SJ) 

N 

Z {N, M) = ^ Jl / im) 5{ni + --- + njv) , (113) 

{n,}»=l 

where f{n) is the one-site weight defined in Eq.®. The sum can be rep- 
resented in the form of the contour integral, 

=/<Si^i|, (U4, 
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where F{z) is the generating function of the one-site weights, 

oo 

F{z) = Y,^V{n). (115) 

Once the partition function is known, it can be used to obtain stationary 
correlation functions. For example the probability P (n) for a site to hold 
n particles is as follows, 

^ / X ... Z(N -l,M-n) 

Another tool characterizing the occupation number distribution is the gen- 
erating function of its moments (e'*'"), which can be also represented in the 
form of contour integral 

(e-) = Y.'-^ (117) 

fe=0 

1 r{F{z)f-^F{ze^)dz 



Z{N,M)J 0-^+1 27ri' 

Another quantity of interest is the probability 'P(n) for n particles to hop si- 
multaneously. The integral representation of the corresponding generating 
function is 

*W-E^»nn)^^/fi^|-, (US, 

where 

oo 

^{x,z) = J2f (") ^" i^P in) + (1 - P(n))) . 

n=0 

This yields for example the following expression for the average flow of 
particles 



1 , / X 1 / \F(z)]^ z dz , , 



The subscript A'' specifies that the expression is valid for an arbitrary finite 
lattice. 



33 



B. The saddle point approximation. 



The integral in Eg. (|1 14(1 can be evaluated in the saddle point approx- 
imation. The location of the saddle point z* is defined by the following 
equation. 

p = z*^-^. (120) 
^ F{z*) ^ ' 

This equation can be treated as an equation of state, which relates the 
density of particles p and the fugacity of a particle z* . It is convenient to 
introduce the Helmholtz free energy 

A{N,M) = -hiZ{N,M) (121) 
~ ~N\tlF{z*) + M\tlz*. 

Then one can define the chemical potential 

dA{N, M) 



dM 

and the analogue of pressure 

dA{N, M) 



dN 



= lnz* (122) 

N 



= lnF(z*). (123) 

M 



The stationary correlation functions mentioned above can also be expressed 
in terms of the values defined. 

- The occupation number probability distribution P{n) for n <^ N 

P{n) = /(n)e"^-'^. (124) 

We should note that when n is of order of N this equation loses the validity 
as one should take into account the shift of the saddle point in the integral 
for Z{N — 1, M — n). In other words, when n in given site is large it 
changes the density of particles in the other sites, which in its turn leads 
to the change of the fugacity of particles. 

-The cumulants {n^}^ of the occupation number n are given by the 
derivatives of log(e'''"), Ea. ljllTII . which leads to a standard relation be- 
tween fluctuations and pressure 



- The average flow of particles <j) is given by 
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The 1 /N correction to the average flow can be given in terms of two first 
cumulants of the occupation number {n^}^, 




(127) 



which in its turn can be reduced to a very simple form 

where one should take into account the equation of state (|120ll . The 1/N 
correction to the average flow has been claimed to be universal for the 
KPZ class, depending only on the parameters of the corresponding KPZ 
equation. Below we use it to reexpress the above parameters fci,fc2 obtained 
from the solution of the dynamical problem in terms of the parameters of 
the stationary state only. 

We should note that the criterion of validity of the saddle point approx- 
imation, i.e. smallness of the second term of the asymptotical expansion 
comparing to the first, yields the following range of parameters: 

« ViV. (129) 

Let us turn to the examples. 

- For the case q = 0, the series F{z) expressed in terms of the parameter 
X ~ p/{l — p) can be summed up to the form 

m = T^ + (130) 
1 + A A — z 

In this case the integral in H114|) can be evaluated explicitly resulting in 

Z (». M) = ,1 + A,- ... ( -f:/ ; -A 

(131) 

Using then the formula (|119|) we arrive at the expression of the particle 
flow (f) given in Ht)5|) . 

- In the case of arbitrary q, \q\ < 1, due to speciflc form of the one site 
weights /(n), the function F{z) has the structure of well known q— series. 



35 



where 

n-l 

{a;q)n='[[{l-aq'') (133) 

is a notation for shifted factorial. The series (|132|l can be summed up 
to the infinite product using the g— binomial theorem |3lj |. which is the 
g— analog of the Taylor series of a power law function, 

Fiz)^{l-p)t^^. (134) 

This formula (|134|l turns out to be extremely useful as it allows one to 
write the equation for the saddle point (|120|l in the following simple form. 

p^y^l^ + y (135) 

k=0 ^ k=0 ^ ' 

Note that if we define 

Xc = z*/\, (136) 

then the r.h.s of Ea. (|135|) coincides with gq^\ (xc) in the form H95|) . Thus, 
the equation for the saddle point is nothing but the conical point equa- 
tion, Eq. (|8(J|I , appeared in the thermodynamic solution of the Bethe ansatz 
equations, while the above parameter Xc is proportional to the fugacity of 
a particle z* . 

Using the equality 

9q,X [Z/X) - 

and Ea. H125|l one can obtain the explicit form of the first two cumulants 
of the occupation number 

= y9;A(=7A), (137) 



inX-^^^.'^^^^]. (138) 



Substituting them into Eas. H12t)ll27|) one can make sure that the expression 
for the particle flow coincides with that obtained from the Bethe ansatz 
solution. 

It is remarkable to note that the variance of the occupation number, 
{■n?) , can be directly related to the parameter fc2. 




27r(n2)^ = J27r-^. (139) 
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At the same time the 1/iV correction to the average flow, Ea. H127|) . is given 
by the product fci/c2 such that 

d^<j) 1 \ 

(140) 

Thus, we have expressed the parameters characterizing the fluctuations of 
the particle current in terms of the parameters of the stationary state, which 
characterize the fluctuations of occupation numbers. As the obtaining of 
the latter does not requires an integrability, we expect this relation to hold 
for the general discrete time ZRP belonging to the KPZ class. 

Another interesting relation can be found. Let us consider the free energy 
per site, a(p, z) as a function of the density p and the fugacity z that 
formally can take on arbitrary complex values 

, , A(M,N) , , 
a[p,z)= lim — = — Ini^(z) + yolnz. 

Then the following relation holds 

where Rq (z) is the density of the Bethe roots in the thermodynamic limit, 
Ea. H69(l . Thus the difference |a(p, zi) — a{p, Z2)\ with the values of 21, Z2 
taken at the contour of the Bethe roots gives the fraction of the roots in 
the segment between ziand 22- We should also note that the equation of 
state, which in terms of a{p, z) looks as follows 

da{p, z)/dz = 0, 

is equivalent to the conical point equation Ea. (|71|> . 



C. The occupation number distribution in the deterministic limit. 

As we noted above, the saddle point equation, Ea. (|135() . for the integral 
that represents the partition function Z{M,N), Ea. H114|) . coincides with 
Ea. H8l)|l appeared in the thermodynamic solution of the Bethe equations. 
Therefore, we can directly use the results of the solution of this equation 
obtained in the section to obtain the occupation number distribution in 
the limit St — > 0. 
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In the domain p < and n <^ N we have 
P{0) = 1-p-0{5t), P{1) = p-0{5t) 



P{n) = 6t 



(1-p) (1-pJ- 



(1 - Pc) P 



(1 - g") (1 - L (1 - p) p. 

For 1 the latter formula yields the exponential decay. 

Pc- P 



(141) 

0{6T'^),n> 2 (142) 



(1 - Pc) 



(l-p)Pc 



When the density approaches the critical point from below, {p^~ p) <C 1, 
the cut-off of this distribution diverges proportionally to {p^ — p)~^.Note 
that strictly below p^ the cut-off is always finite, being independent of 5t, 
while P{n) itself is of order of St. 

- In the domain p> p^ and n <^ N we have 



P(0) = 1 - - 0(<5t1/2), P(l) = - 0{5t''^) 
and for n > 2 



P(n) = 5t 



(1 - 9") (1 -«"-!) 



1 



(1 - Pc) (P - Pc) 



(143) 

■0((5t3/2). 

(144) 



For large n the formula for P(n) also yields the exponential decay 

5t 



Pin) 



(1 - Pc) 



exp 



St 



(1 - Pc) (P - Pc) 



(145) 



However, in this case the cut-off depends on 6t diverging as St~^^'^ as St 
goes to zero. 

- When p = Pc and n<^ N we have 



P(0) = 1 - - 0{STy% P(l) = p, - OiST^') 



(146) 



and for n > 2 



P(n) 



'^T(l-pJ-^ 
(l-9")(l-g"-i) 



(5r 



Pc(l - Pc)^ 



1/3- 



+ (jr^/^) . (147) 



In this case the cut-off of P(n) at large n is of order of St 



P(n) 



St 
(1 - Pc) 



exp 



Pc(l - Pc)^ 



1/3- 



(148) 
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One can sec from these formulas that below the critical point, the density 
of MP sites at the lattice vanishes proportionally to St 

M r 2 

while their average occupation number {n) ^jp is finite, 

both increasing as {p^ — p)~^ as p approaches p^. Above p^ the density of 
MP sites vanishes as V&r, 



ISt {p ~ p^) 

PMP{P> Pc) - \l ) ' (l^l) 

which is much slower than in the subcritical region. Their average occupa- 
tion number diverges as 5t~^/'^ 



(1-Pc)(p-Pc) 
5t 



{n),,p ip > pj ^/^ '^-^r ■ (152) 



Exactly at the critical point we have 



^^2/3^1/3 

PMpiP = Pc) - nTJ: (153) 

(1 - Pc) 



IMP 



{P = pj^(^)'^'(l-pj^/^ (154) 



One can see that the fraction of particles in MP sites, i.e. Pmp {^) mp^ 
above the critical point is exactly equal to {p — p^), which preserves the 
density of particles in MP sites equal to p^. 

The above formulas are derived under the suggestion n <ti N. For n of 
order of N one should take into account the shift of the saddle point in 
the numerator of Eq. I|116l) . The criterion of the applicability of the saddle 
point method (|129|l suggests 

St^/'^N > 1 (155) 

for p > p^ and 



St^^^N > 1 



(156) 
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FIG. 3: The probability distributions P{n) of the occupation number n for two 
values of density and q = 0.5, 5t = 10"^, iV = 100 

for p — p^. Under these conditions the cutoffs of the distributions 1146t 
I148|) are much smaller than iV, and the probability of MP sites with n ^ N 
is exponentially small in N. One can address the question what happens 
beyond the domains H155ll5t)|) . The extreme limiting case corresponds to 
the situation of the ASAP on a finite lattice, when St is so small that one 
can neglect any quantities of order of 5t^ . In this case there is a tendency of 
creation of one MP site with the occupation number N(p— p^) and the life 
time that grows exponentially with the system size N [35j . The situation in 
the intermediate region should be addressed separately. Direct calculation 
of the occupation number distribution shows a peak in the middle, which 
grows as St approaches zero (see Fig (O). Analysis of the finite lattice 
behaviour is beyond the goals of present article. 



We have presented the Bethe ansatz solution of the discrete time ZRP 
and ASEP with fully parallel update. We found the eigenfunctions of the 
equation for configuration-dependent generating function of the distance 
travelled by particles. The eigenfunctions of this equations were looked for 
in the form of the Bethe function weighted with the weights of stationary 
configurations. We started with the ZRP with arbitrary probabilities of 
the hopping of a particle out of a site, which depend on the occupation 
number of the site of departure. We inquired which restriction on the 
probabilities are imposed by the condition of the Bethe ansatz solvability. 
As a result we obtained the two-parametric family of probabilities, the two 
parameters being p, the hopping probability of a single particle, and q, the 
parameter responsible for the dependence of probability on the occupation 
number of a site. The model turns out to be very general including as 



VII. SUMMARY AND DISCUSSION. 
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particular cases g-boson asymmetric diffusion model, phase model, drop- 
push model and the ASAP. By simple arguments the ZRP can be also 
related to the ASEP-like system obeying exclusion interaction. In such 
system the hopping probabilities depend on the length of the queue of 
particles next to a hopping particle. The particular case of the model is 
the Nagel-Schreckenberg traffic model with f,nax = 1- We obtained the 
Bethe equations for both cases. 

The further analysis was devoted to the calculation of generating function 
of cumulants of the distance travelled by particles. In the long time limit it 
is given by the maximal eigenvalue of the equation discussed. We obtained 
the largest eigenvalue in all the phase space of the model. First, the case q = 
corresponding to the Nagel-Schreckenberg traffic model was considered, 
which due to special factorized form of the Bethe equations can be treated 
exactly for an arbitrary lattice size. The result is consistent with that 
for the continuous-time ASEP and the stationary solution of the Nagel- 
Schreckenberg traffic model. For the general q we obtained the maximal 
eigenvalue in the scaling limit. It has the universal form specific for the 
KPZ universality class. Again, the model dependent constants obtained 
reproduce correctly all the known particular cases. 

We carried out detailed analysis of the limiting case p ^ 1. The phase 
transition which takes place in this limit was studied. We shown that below 
the critical density, the flow of particles consists of almost deterministic 
synchronous jumps of all particles. Above the critical density the new 
phase appears. The finite fraction of all particles gets stuck immobile at 
the vanishing fraction of sites. The fluctuations of the particle current 
become singular, nonanalytic in (1 — p). 

In terms of the associated ASEP the transition studied is the analytic 
continuation of well-known jamming transition in the Nagel-Schreckenberg 
model with v,„ax — 1, which corresponds to a particular case of our model. 
The value of the critical density in that case, = 1/2 (or — 1 for 
associated ZRP), follows from the particle hole symmetry. Switching on 
the interaction between particles breaks this symmetry and as a result 
decreases p^. 

From the point of view of the hydrodynamics of the particle flow the 
reason of the jamming is the vanishing of the velocity of kinematic waves 
Vkin = d(f)/dp, which are responsible for the relaxation of inhomogeneities. 
As was argued in 37], the situation Vkin = leads to the appearance of 
shocks in the ASEP interpretation or the MP sites with large occupation 
numbers in ZRP. When p > p^ and p ^ 1, Vkin behaves as follows 



\/i -p p: 




1 1/2 



d<f)/dp 



2 



1 /3 

while in the close vicinity of the critical point, \p — ^ {1 — p) , the 
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1/3 

velocity of kinematic waves vanishes as (1 — p) . Similar transition also 
was considered in the bus route model |23 ■ In that case however the rates 
u{n) of hopping of a particle out of the site occupied by n particles depend 
on the external parameter A, such that the limits A — )■ and n ^ oo are 
not commutative. In our case the transition owes to taking a limit in only 
one hopping probability p(l) and to the discrete parallel update. 

We also obtained the canonical partition function of the stationary state 
of the discrete time ZRP. It is expressed in terms of the q— exponential 
series. We apply so called g— binomial theorem to evaluate the partition 
function in the saddle point approximation. Using the partition function 
formalism we obtained the formulas for physical quantities characterizing 
the stationary state of the model such as the occupation number distri- 
bution and the average flow of particles, the latter confirming the Bethe 
ansatz result. We observed several curious facts, which reveal intimate 
relation between the saddle point approximation applied for the partition 
function formalism and the thermodynamic limit of the Bethe ansatz solu- 
tion. We found that two constants which define the behaviour of the large 
deviation function of the particle current can be expressed in terms of the 
cumulants of the occupation number of a site and the fugacity of particles. 
Furthermore, the analogue of the Hclmholtz free energy considered as a 
function of arbitrary complex fugacity plays the role of a "counter" of 
the Bethe roots. The equation of state, which relates the density and the 
fugacity, coincides with the conical point equation, H71|) .which fixes the po- 
sition of the Bethe roots contour. All the correspondences found seem not 
accidental. It was claimed in [s^l that the large deviation function for the 
particle current should be universal for all models belonging to the KPZ 
universality class. On the other hand, once the density of the Bethe roots is 
obtained, the derivation of the nonuniversal parameters is straightforward. 
It would be natural to expect that the derivative of the free energy would 
serve as a generalization of the density of Bethe roots for the cases where 
the Bethe ansatz is unapplicable. Currently we do not have any explana- 
tion of the relations found. However, if they do exist, clarifying of their 
internal structure could give a way for study of a variety of systems which 
do not posses the Bethe ansatz integrability. We leave these questions for 
further work. 

Another possibility of continuation of current study is as follows. In 
the particular case, g = 0, of our model, which corresponds to the Nagel- 
Schreckenberg traffic model with Wmax = 1, the Bethe equations have spe- 
cial factorized form. Usually, the consequence of such factorizability is the 
possibility of time-dependent correlation functions in the determinant form. 
This program was realized before for the continuous time ASEP and 
the ASEP with backward ordered update 0| . We expect that the parallel 
update case is also treatable. 
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